sor_solver <- function(H, q, lb, ub, u0 = (lb + ub) / 2, max.steps = 500) {
  n <- ncol(H)
  u1 <- u0
  Hui <- H * matrix(u1, n, n, byrow = T)
  Hu <- H %*% u1
  for (j in 1:max.steps) {
    numerator <- q - Hu
    for (i in 1:n) {
      u1[i] <- u0[i] + (q[i] - Hu[i])/H[i, i]
      u1[i] <- max(min(u1[i], ub[i]), lb[i])
      Hui_ <- H[, i] * u1[i]
      Hu <- Hu - Hui[, i] + Hui_
      Hui[, i] <- Hui_
    }
    if (norm(u1 - u0, type = "2") < 0) {
      break
    } else {
      u0 <- u1
    }
  } 
  return(list("x" = u1))
}

H <- manysvms:::calculate_svm_H(Xe %*% t(Xe), y)
q <- matrix(1, nrow(X), 1)
lb <- matrix(0, nrow(X), 1)
ub <- matrix(1, nrow(X), 1)


u <- sor_solver(H, q, lb, ub, lb, 50)$x
t(Xe) %*%  (u*y)

system.time(model <- pin_svm2(X, y, max.steps = 20000, tau = 1))
plot(model)
coef(model)

system.time(model2 <- manysvms::pin_svm(X, y, max.steps = 320000, eps = 0, tau = 1))
coef(model2)
plot(model2)

system.time(model <- hinge_svm2(X, y, max.steps = 40))
plot(model)
coef(model)

system.time(model <- hinge_svm(X, y, max.steps = 800))
plot(model)
coef(model)

